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Abstract. - We study the response of a granular system at rest to an instantaneous input of 
energy in a localised region. We present scaling arguments that show that, in d dimensions, the 
radius of the resulting disturbance increases with time t as t™, and the energy decreases as t~ ad , 
where the exponent a = l/(d + 1) is independent of the coefficient of restitution. We support 
our arguments with an exact calculation in one dimension and event driven molecular dynamic 
simulations of hard sphere particles in two and three dimensions. 



Granular systems, predominantly characterized by dis- 
sipative collisional dynamics, are ubiquitous in nature and 
exhibit a wide variety of very rich and striking physical 
phenomena [1]. Although many experimental studies have 
captured the complexity of these systems by studying phe- 
nomena ranging from clustering instability, co-existence of 
phases to non-Maxwellian velocity distributions (see [1,2] 
for reviews), the theoretical understanding of these sys- 
tems is far from complete (see [3-5] for reviews). Hence, 
it is imperative to study simple models that capture some 
distinctive features of the system, yet are amenable to 
analysis. 

A model that has attracted considerable attention in 
, the past is the freely cooling granular gas, where the par- 
ticles move ballistically and lose energy only through in- 
elastic collisions [6-17]. Starting from a homogeneous spa- 
tial distribution of particles with velocities drawn from 
a normalizable distribution function, simulation studies 
show that after an initial regime when energy decays as 
E t ~ t~ 2 (Haffs law) [18], clustering instability sets in [7]. t 1 
The long time behavior of the system is universal: the 
energy decays algebraically with an exponent which de- 
pends on the dimension but not upon the coefficient of 
restitution [11-13,17]. The exponent is known analyti- 
cally in one dimension through a mapping to the Burgers 
equation (E t ~ t~ 2 / 3 ) [10,19]. In higher dimensions, the 
exponents obtained from the analogy to Burgers equation 
{E t ~ t~ d/2 ,d > 2) [12] differ from that obtained from 
mean field scaling arguments (E t ~ ^-2d/d+2-j jgj an( j f rom 
simulations of the Boltzmann equation [13,20], leading to 
an uncertainty in the precise value of the exponents in two 



and higher dimensions. 

In this paper, we consider a simple and tractable model 
of a cooling granular gas where the particles are initially 
at rest and the system is perturbed by imparting momen- 
tum to a single particle. This in turn leads to motion 
of other particles by inter-particle inelastic collisions, and 
the particles cluster to form a nearly spherical shell that 
propagates radially outwards in time [see Fig. Q] (a)]. Us- 
ing scaling arguments and numerical simulations, we show 
that the scaling behaviour of energy and the radius of the 
disturbance with time is independent of the coefficient of 
restitution. The results obtained from scaling arguments 
are confirmed by an exact calculation in one dimension 
and event driven molecular dynamics simulations in two 
and three dimensions. 

The corresponding problem when collisions are elastic is 
the classic Taylor-von Neumann-Sedov problem of shock 
propagation following a localized intense explosion [21]. In 
this case, the particles remain homogeneously distributed 
see Fig.[T](b)] and the exponents can be obtained by sim- 
ple dimensional analysis [22], while the scaling functions 
can be calculated exactly following a more detailed anal- 
ysis [21,23]. The simulations and scaling arguments for 
a hard sphere model with elastic collisions were recently 
done in Ref. [24]. Signal propagation has also been stud- 
ied in excited dilute granular gas [25] as well as in dense 
static granular material (see [26] and references within). 

Our model consists of a collection of monodisperse hard 
spheres (in simulation we have taken 2.5 x 10 5 and 2 x 10 6 
particles in two and three dimensions respectively) of fi- 
nite diameter (unity in simulation) distributed randomly 
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Fig. 1: (Color online) Shown are the positions of particles that 
have undergone at least one collision, following input of energy 
at (500, 500) for (a) the inelastic case (r = 0.1) at times t — 
5000, 10000, 20000, 50000, and (b) elastic case (r = 1.0) at time 
t = 25000. 



in space such that no two particles overlap (in simulation 
the number density n = 0.25 in both two and three dimen- 
sions). Periodic boundary conditions are implemented in 
all directions. All the particles are initially at rest. A 
single particle is chosen at random and given a velocity 
of unit magnitude along a random direction. The particle 
motion is ballistic till it collides with other particles. The 
collisions conserve momentum and the velocities change 
deterministically according to the following collision rules: 
if the velocities before and after collision are Ui, 112, and 
vi, V2 respectively, then 



a is a scaling exponent. Then, 

±a — l 



Vi 



ltd 



ad+2a-2 



(2) 
(3) 
(4) 



The above relations hold good for both elastic and in- 
elastic collisions. We now analyze the two cases separately. 
For the elastic gas, energy is a constant of motion. This 
implies 

a = dh> r=1 > (5) 

coinciding with the results for one and two dimensions in 
Ref. [24] and for three dimensions in Ref. [22] . 

For the inelastic case, there is one unknown exponent a 
which is determined by the following argument. A short 
time after the initial perturbation, the particles that have 
undergone at least one collision concentrate themselves 
into a narrow band. Though the data shown in Fig. Ufa) 
is for r — 0.1, clustering is seen for all r < 1. Due to 
this spatial structure, there is no radial momentum trans- 
ferred from particles at a certain angle to those that are 
diametrically opposite, or in other words, the radial mo- 
mentum is conserved. The radial momentum carried by 
the particles in a small solid angle dfl scales as vtRfdil. 
The conservation law implies that vtRf ~ const, or equiv- 



alcntly, vt ~ Rt 
immediately obtain 



-ad 



Comparing with Eq. (|2|), we 



Vi, 2 = U1.2 - e[n.(ui j2 - u 2 ,i)]n, 



(1) 



where r — 2e — l(0<r<l)is the coefficient of restitution 
and n is the unit vector directed from center of particle 1 to 
center of particle 2. Thus, the tangential component of the 
relative velocity remains unchanged, while the magnitude 
of the longitudinal component is reduced by a factor r. 

For r < 1, the system undergoes inelastic collapse in 
which infinite collisions take place in finite time [27] . This 
computational difficulty is avoided by making the colli- 
sions elastic when the longitudinal relative velocity is less 
than a cutoff velocity 5 [11]. This qualitatively captures 
the experimental situation where r is seen to be a function 
of the relative velocity [28,29]. In our simulation, we set 
5 = 10~ 4 . 

Consider now the result of a typical simulation [see 
Fig. [Ha)]. Let Rt be the typical radius of the shock pro- 
file, v t the typical speed, N t the number of active particles 
(particles that have undergone collisions), and E t the to- 
tal kinetic energy at time t. These quantities are related 
to each other through simple scaling relations. The speed 
i>t is related to Rt as Vt ~ dRt/dt. The number of parti- 
cles that have undergone collisions is proportional to the 
volume swept out by the disturbance: Nt ~ Rf, where d 
is the dimension. Energy is then given by Et ~ Ntvf. 

We look for scaling solutions of the kind R t ~ t a , where 



1 



< 1. 



(6) 



In one dimension, the above scaling result can be 
checked by a simple calculation. Consider the sticky limit 
r = 0, when the particles coalesce on collision. Let par- 
ticles of unit mass be initially placed on a lattice with 
spacing a. Let the particle at the origin be given a ve- 
locity vo to the right. When this particle collides with its 
neighbor, it coalesces with it. The mass of this compos- 
ite particle after m collisions is then m, and its velocity, 
given by momentum conservation, is v m = vo/m towards 
the right. The time taken for m collisions is given by 



m— 1 

v - a 



1 = 

am(m — 1) 



2vn 



(7) 
(8) 



At large times, m w y // 2vot/a. But m is identical to Nt 
and Rt, which by definition scales as t a . This gives a = 
1/2, consistent with that obtained by setting d = 1 in 
Eq. ©. 

In two and three dimensions, the scaling arguments are 
tested numerically using event driven molecular dynamics 
simulations [30] . The data presented is averaged typically 
over 100 different initial realizations of the particles. All 
lengths are measured in units of the particle diameter, 
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Fig. 2: The anisotropy index A(t) of the band in two dimen- 
sions is plotted as a function of time t for different values of 
the coefficient of restitution r. A(t) converges to a value less 
than one for all r. 



and time in units of initial mean collision time 1/ '(von 1 ^), 
where vq is unity in the simulations. We hrst check the va- 
lidity of the assumption of a single length scale Rt- Fig. [2] 
shows the variation in two dimensions of the anisotropy in- 
dex A(t) with time, where the anisotropy index is given by 
A(t) = ([(Ai-A 2 )/(Ai+A 2 )] 2 ), Ai, A 2 being the eigenvalues 
of the moment of inertia tensor [31]. If the transverse and 
longitudinal radii scale differently with time, then A(t) 
should converge to unity at large times. However, A(t) is 
found to converge to a constant less than one for all r. For 
r = 1, Ait) converges to zero at large times. We conclude 
that though the shape of the front is anisotropic for r < 1 , 
all length scales scale identically with time. 

We check the scaling relations Eqs. ([3]), ((3]), and © by 
measuring the mean number of active particles (Nt) and 
the mean total kinetic energy (E t ) as a function of time. 
In two dimensions, the scaling argument gives (N t ) ~ t 2 ^ 3 , 
(E t ) ~ t~ 2 / 3 , while in three dimensions (N t ) ~ i 3 ^ 4 , 
(Et) ~ i -3 / 4 . In Fig. [3]^ a) and (b), we show the variation 
with time of Nt and Et in two and three dimensions. For 
larger r, it takes longer time to reach the scaling regime. 
This crossover time $p reflects the transition of the par- 
ticles from the initial homogeneous spatial distribution to 

tc diverges in the elastic 
(1 — r 2 )^^ 1 , where cf>\ w 2.25 in two dimen- 



the clustered state. We hnd that t c 
limit as tP - n 
sions and <j>\ w 3.0 in three dimensions. In addition, at 
large times, the system crosses over to the elastic regime 
when v t ~ 8. This crossover time scales as ~ 5~^ 2 
where 4> 2 = 1/(1 - a) [3/2 in d = 2 and 4/3 in d = 3]. 
Within these limitations, the numerical data shows good 
agreement with the theoretical prediction shown with solid 
lines. 

We check the scaling relations for R t and vt by studying 
the radial and the velocity distribution function. The ra- 
dial distribution function P(R, t) measures the mean num- 
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Fig. 3: Simulation results for the (a) the mean number of ac- 
tive particles (Nt) and (b) the mean kinetic energy (Et) as a 
function of time t. In both the plots, the top three curves cor- 
respond to three dimensions and the bottom three curves cor- 
respond to two dimensions. The different data correspond to 
the coefficients of restitution r = 0.1(O), 0.5(A), 0.8(D). The 
solid lines have exponents obtained from scaling theory. The 
data have been shifted for the sake of clarity. 



ber of active particles at a distance R from the center of 
mass of the active particles at time t. The velocity dis- 
tribution function P(v,t) measures the probability that 
a randomly chosen active particle has speed v at time t. 
These distribution functions should be a function of a sin- 
gle scaling variable: 



P(R,t) = t- a h(Rt- a ), 
P(v,t) = ^-"Uvt 1 -"), 



(9) 
(10) 



where fx and f% are scaling functions. These scaling 
collapses are verified numerically in two dimensions (see 
Fig. and in three dimensions (see Fig. 0. The data 
shown is for one value of the coefficient of restitution 
(r = 0.1), but the same is observed for other values of 
r. The scaling function f 2 (vt 1 ~ a ) decays exponentially at 
large speeds v. Such non-Maxwellian behaviour is typi- 
cal of granular systems [32-34] . We also observe that the 
faster particles are in the inside edge of the collapsed band, 
thus making the bands stable. 

We also studied the structure of the collapsed bands. 
For that, the packing fraction of the particles in the bands 
was numerically calculated by dividing the space into cells 
of linear length 10, and counting the number of particles 
in each cell. For all r < 1, the typical packing fraction seen 
at large times ranges from 0.78 — 0.82 in two dimensions. 
This value is very close to 0.84, the packing fraction of ran- 
dom close packed structures seen in jamming of frictionless 
spherical particles [35]. For r — 1, the packing fraction is 
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Fig. 4: Results in two dimensions for (a) the radial distribu- 
tion function P(R, t) and (b) the velocity distribution function 
P(v,t), when scaled as in Eqs. © and (|10[) with scaling ex- 
ponent a — 1/3. The scaling collapse has been obtained for 
times t = 25000(0), 37500(A), and 50000(D). The coefficient 
of restitution is r = 0.1 



Fig. 5: Results in three dimensions for (a) the radial distribu- 
tion function P(R, t) and (b) the velocity distribution function 
P(v,t), when scaled as in Eqs. <(9j and (|10|l with scaling ex- 
ponent a = 1/4. The scaling collapse has been obtained for 
times t = 30000(0), 50000(A) and 75000(D). The coefficient 
of restitution is r = 0.1. 



~ 0.47, showing that the particles are very loosely packed. 

To conclude, we studied the problem of shock propaga- 
tion in granular (inelastic) systems and obtained scaling 
solutions for the problem. In one dimension, the exact 
result for the sticky limit (r = 0) corroborated the scal- 
ing solution. In two and three dimensions, we verified 
our results using event driven molecular dynamics simu- 
lations. Our analysis showed conclusively the universality 
(r-independence) of the scaling solutions and its depen- 
dence only on the spatial dimension. We retrieved the 
earlier results for the classic Taylor-von Neumann-Sedov 
problem corresponding to the elastic limit (r = 1). For 
r < 1, we obtained an explicit expression for the scaling 
exponent in the late time cooling which has hitherto re- 
mained inconclusive for the related problem of the freely 
cooling granular gas. 

The model discussed in this paper also has experimental 
significance. Direct experiments on freely cooling gas are 
difficult due to friction and boundary effects. Recent ex- 
periments reproduced the energy decay law in the homoge- 
neous cooling regime [36], but not in the clustered regime. 
The boundary effects will be eliminated if the granular gas 
is initially at rest, making the problem discussed in this 
paper more easily reproducible in the laboratory. 
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